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1. Introduction 


Geothermal energy is considered a strategic resource in many 
countries, even if its use appears to be often marginal in the 
national energy systems. Its continuous operating mode distin- 
guishes geothermal energy from the other renewable sources, 
intermittent or stochastic. Majority of the geothermal sources 
worldwide are medium-low enthalpy type (water dominant, at 
temperature lower than 150°C and pressure below 15 bars). 
Stefansson, |1] estimated that more than 70% of the geothermal 
resources available in the world are water dominated fields, at 
temperatures under 150°C. The binary cycle technology with 
Organic Rankine Cycle (ORC) appears to be the most efficient 
and convenient solution for such a kind of resource [2]. Binary 
power plants are now objects of wide attention by energy markets, 
although their diffusion is still made difficult by a lacking technol- 
ogy standardization and due to the quite high specific costs [3,4]. 
The great variability of the resource characteristics worldwide is 
one of the possible reasons. The proper matching between the 
reservoir capability and the plants parameters (power size, extrac- 
tion/reinjection rate) is a critical key point. In power plants using 
dry-steam (high enthalpy) geothermal resources, pressure and 
temperature reduction can be compensated by an increase of the 
mass flow rate. In case of binary plants, a variation of the resource 
properties (Tzeo, Pgeo) could also lead to a fast end of life of the 
plant. The first and most important activity to design a geothermal 
energy plant is an accurate investigation of the geothermal 
potential assessment, as well as the prediction of reservoir 
response at given industrial exploitation configurations. For these 
reasons, a multidisciplinary approach to the problem of exploita- 
tion of geothermal fields (in particular at medium-low tempera- 
ture) is necessary. Thermal engineering, geochemistry, geophysics, 
and reservoir engineering are the fields involved in this technique 
(Fig. 1). The authors have diffusely discussed this topic in a recent 
paper [5]. 

Numerical simulation is a fundamental and strongly interacting 
instrument for plant design [6]. Different approaches are here 
considered with reference to several case studies of geothermal fields, 
which are reviewed and discussed. The perspectives of numerical 
simulation of geothermal reservoirs as support to the design and 
sizing of geothermal plants are also outlined. Simulation can be very 
important in order to define and progressively modify the manage- 
ment strategy of the geothermal field. Construction of the numerical 
model must be supported by a detailed knowledge of the spatial 
distribution of the properties of the reservoir: the accuracy in the 
definition of the dataset is fundamental for the construction of the 
model. The model is then enriched by including the database of 
historical data collected during exploration. 

The results obtained depend a lot on the accuracy level of the 
input data. The model will be much more accurate if as much 
details as possible are known about the geological properties of 
the rocks (effective porosity, density, specific heat, permeability, 
thermal conductivity), thermophysical properties of the fluid 
(specific heat, density, thermal conductivity), fractures pattern 


Thermodynamics 
Energy engineering 


Geochemistry Reservoir 
engineering 


Geophysics 
Fig. 1. The multidisciplinary approach proposed, with the connections between the 
three areas involved. 


and layout, natural recharge of fluid, geothermal boundary 
conditions. 

The numerical model of a geothermal reservoir is very impor- 
tant both for the definition of the geothermal potential assessment 
and of the reinjection strategy. The geothermal potential of a 
particular area means the definition of temperature (Tgeo) and 
pressure (Pgeo) Of the geothermal fluid and also of the maximum 
mass flow rate (Mge) that can be extracted maintaining the 
thermal properties of the reservoir and of the geofluid constant 
for a long time. Concerning the reinjection strategy, it is necessary 
to take into account the circulation model of the fluid in the 
regional area considered [7]. A general methodology for the 
reinjection technologies is not properly available in literature, 
the optimal strategy is in fact site-dependant, as the potential 
assessment itself. Interesting discussions on this particular topic 
are reported by Sigurdsson et al. [8], and recently by Kaya et al. [9]. 

The main task of potential assessment and sustainable plants 
design is the optimization and enhancement of the resource 
durability (Vaccaro [6]). Interesting contribution on the definition 
and evaluation of the sustainability and renewability of the 
geothermal energy uses are available in the recent paper by 
Hahnlein et al. [10], together with the paper by Axelsson [11]. 
Particularly in case of innovative geothermal utilizations (like for 
example EGS) the long-term consequences on the environment 
are not completely known yet. The same argument can be then 
referred to the renewability of the resource itself, which is directly 
dependent on the type and rate of utilization. The renewability 
(and sustainability) reference level can vary, as one can adjust the 
energy system size and extraction rate according to an acceptable 
durability level [11]. Also in case of direct heat utilization (e.g. 
district heating) some strategy remarks should be pointed out. In 
the recent paper by Fox et al. [12] the renewable capacity of deep 
systems is assessed and discussed in order to elaborate a rotating 
utilization strategy. 

The numerical simulation of geothermal reservoirs is a well 
known topic and has already been an object of investigations and 
reviews (e.g. O'Sullivan [13]). Unfortunately, till now the use of 
numerical simulation has not faced any direct connection with the 
energy systems analysis. A proper prediction should deal with the 
changes of the different parameters in response to given mass flow 
rate extraction and reinjection (corresponding to the specific 
energy strategy). It is evident that a key role is assigned here to 
the numerical simulation of the reservoir, as compared to other 
reservoir engineering aspects (wells siting, fluid losses, tracer test). 
This ambitious task seems to be not of specific interest in most of 
the analyses carried out in the past, so the authors would like to 
review the recent developments in the field of numerical simula- 
tion of geothermal fields, focusing their attention on the particular 
use of such an important instrument for the sustainable design of 
geothermal plants. 


2. Numerical simulation of geothermal reservoirs: strategic 
role for the design of geothermal plants 


The numerical simulation of a geothermal reservoir is a well 
known field of research in the literature and it has already been an 
object of accurate review analysis and methodological overview 
[13-17]. 

Two main goals can be identified: history matching and 
forecast of future scenarios (consequent to the exploitation of 
the reservoir). History matching is usually done to check the 
reliability of a model and evaluate the sustainability level in 
retrospect. It is the analysis of an exploitation history according 
to the data log until present time or during a particular time 
interval. This also allows to check the numerical model in a 
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“feedback loop” and calibrate or adapt it to what is reported in the 
history data (in order to improve future scenarios forecast). This is 
practically done by assigning temperature and mass flow rate of 
both the geofluid extracted and reinjected, and any other useful 
historical data recorded that can be translated into thermophysical 
parameters or boundary/initial conditions. Some phenomena 
which affect the behavior of field utilization could not be put into 
a simulation in a proper way, without knowing about them from 
the history (natural recharge, natural change of the pathways of 
circulation into the rock formations, losses of pressure, etc.). If a 
model is considered to be reliable, it can be used to study how the 
durability of the resource changes depending on different extrac- 
tion rates, reinjection temperatures, wells siting, fractures (also 
induced). A productive strategy can be based on the results of a 
model simulation. This is true for both power plants and thermal 
energy direct uses. Some general and macroscopic aspects of the 
application of numerical simulation to the geothermal energy 
study are listed here: 


è equations describing the phenomena considered (circulation, 
energy/mass transport); 

estimation of the main thermophysical parameters; 

boundary conditions (BC); 

geothermal potential assessment; 

coupling between the power/thermal utilization and reservoir. 


2.1. Characterization of the geothermal source potential 


Potential assessment is a fundamental step of a geothermal 
project and its final goal is the sustainable utilization of the 
resource. It involves the complete characterization of the field, 
energy stored, maximum fluid rate, useful temperatures and 
chemical parameters of the fluid (to determine the minimum 
temperature for reinjection). This evaluation is surely important 
for each kind of water dominant geothermal field, but mainly in 
case of moderate temperature geothermal fields in which the 
installation of an ORC plant is programmed. During the running of 
a plant the productivity (and wells deliverability) can show some 
remarkable variations (in terms of flow rate, fluid chemistry and 
specific enthalpy of geofluid). These changes can be addressed to 
reservoir pressure decline because of excessive fluid production. 
The reservoir fluid pressure is also another important parameter, 
linked both to the productivity of the well and the scaling 
phenomena. High pressure in the geofluid pipes can keep the 
scaling phenomena under controlled rates. 


2.2. Reinjection strategy and geofluid chemistry 


Only if reinjection is practiced one can say that geothermal 
energy is being used as a renewable energy source. The practice of 
reinjection avoids temperature and pressure decline in a geother- 
mal field. For the binary power plants it is a basic approach for 
resource management and it appears to be compulsory. The task is 
the optimization and enhancement of the durability of the 
resource: the argument has been diffusely analyzed in [6]. The 
objectives that this practice has to achieve are essentially the 
efficient restitution of the geofluid to the reservoir in order to 
optimize the recharge in terms of enthalpy and flow rate, and 
choosing the correct siting and depth of the wells to guarantee a 
sustainable use of the resource. The optimum reinjection strategy 
is a quite complex task and it strongly depends on the type of 
geothermal system (Kaya et al. [9]). In general the minimum 
temperature values are in the range between 60 °C and 80 °C. 

The siting of production and reinjection well and their mutual 
interference are the main issues of the strategy. To give a trivial 


example: if they are too far, the recharge could occur in a long time 
interval (longer than the plant lifetime). If they are too close a cold 
fluid short-circuiting could occur. A correct reinjection strategy has 
to take into account the thermodynamic and chemical problems of 
the scaling phenomena, which in many cases increase with low- 
ering of the temperature. Early study of the geothermal system is 
fundamental, in order to avoid the worst consequences of fouling, 
corrosion and clogging of the parts of the plant, of the pipings, and 
the “tapping” of the wells. 


3. Numerical simulation of geothermal reservoirs in the 
integrated approach to sustainable design 


Numerical simulation of geothermal reservoirs allows us to 
understand the hydrogeological behavior and heat transport into 
the reservoir under a certain utilization rate. It is possible to study 
the geothermal reservoir by solving the balance equations of mass, 
momentum and energy in the particular volume in which hydro- 
thermal circulation of fluid occurs. The hydrogeological issues 
linked to the geothermal exploitation (knowledge of the geological 
structures and of the groundwater system) must be connected 
with the engineering tasks of the design and optimization of the 
energy conversion system. The constitutive laws are peculiar for 
each kind of reservoir, while the numerical analysis issues are also 
important in order to achieve a reliable solution. The development 
of the numerical model itself has to follow two main directions: 


(1) The unperturbed (or undisturbed) natural state. 
(2) The utilization scenarios (during the exploitation). 


Different phases in the development of the model can be 
identified. A first “block-structure” has to be built together with 
the dataset of the parameters which best fit what it is expected by 
the conceptual model. This first step model should reproduce: 


Collecting and interpretation 
of the most meaningful data 


Conceptual model 


Numerical model: geometry, 
equations, BC 


Iterative calibration 
and refinement of 
parameters and BC 
that best fit the 
conceptual model 


Modelling of the unperturbed 
state to verify the model 
stability 


Simple models 


Simulation of the different 
production scenarios 


3D models (more 
complicated) 


Fig. 2. Conceptual scheme for the elaboration of a numerical model of a geother- 
mal reservoir. 
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© Geological structure of the reservoir. 

© Geometrical features of wells and fractures. 

e Hydraulic, Thermal, Mechanic and Chemical conditions 
(HTMC). 


3.1. Conceptual model and numerical model 


The conceptual scheme for the realization and simulation of a 
model of a geothermal reservoir is represented in Fig. 2 (modified 
after [14]). If a numerical model is realized in a multidisciplinary 
team or environment it happens that who is in charge of building 
up the model is not totally aware of the conceptual model of the 
field. It then becomes necessarily a work of a team in which 
different geothermal backgrounds and experiences are shared. An 
important step of “interpretation” and appropriate data collection 
must be pursued firstly, for the elaboration of the conceptual 
model (Huenges [15]). This first step model should reproduce: 
geological structure of the reservoir; geometrical features of wells 
and fractures; Hydraluic, Thermal, Mechanic and Chemical condi- 
tions (HTMC). 

The model should then pass a further step of calibration and 
refinement. It is an iterative process, in which the parameters and 
boundary conditions should be adjusted according to the con- 
ceptual and physical model previously elaborated and to the 
uncertainty level of part of the database used. Some thermophy- 
sical properties (permeability, thermal conductivity, porosity, 
specific heat capacity) change with siting, depth and hydrothermal 
alteration. Their value can be adjusted and the mesh can be refined 
at this step. 

A first model of the unperturbed state is the result of the 
attempts described here. It should be run for a long simulation 
time interval (10°-10° years) in order to verify the model stability 
and convergence. In case of strong uncertainty about both the heat 
transport phenomena and geophysical data, simple 2D models (or 
lumped parameters models) could be firstly run. Exploitation and 
energy utilization scenarios can be then run starting from the 
unperturbed state simulation as initial conditions. The renew- 
ability assessment and durability of the resource have to be results 


of the scenarios simulation. In particular, temperature and pres- 
sure should be kept stable in the reservoir as much as possible. If 
chemical properties and saturation curve of the specific geofluid 
mixture are known, scaling and chemical deposition phenomena 
can be also introduced in the calculations, in a multi-physics 
simulation environment. The models can couple different trans- 
port equations, referring to mass, heat and chemicals. It is 
important to remark that in the results of the simulation some 
effectively useful data for the plant design must be extracted. It is 
evident that some of the geophysical and general results are not 
directly necessary or relevant for who is in charge of designing the 
plant. A close interaction between who elaborates the numerical 
model and who designs the utilization plant would be needed 
(according also to the considerations and remarks about inter- 
disciplinary work and sustainability assessment discussed in the 
previous chapters). 

From a practical point of view, building a numerical model of a 
geothermal reservoir is not that easy. With respect to other 
engineering or science systems, a reservoir cannot be practically 
measured in all its features. Its evolution can be studied, and a 
model can represent a way of behavior prediction only if it is based 
on reliable data. In Fig. 3 a conceptual scheme about how the 
development of numerical models can help the sustainability 
assessment is shown. The size of the geothermal area must be 
known, the geometric domain must be big enough to comprehend 
every interesting mass/energy transfer boundary or spot, but also 
small enough to reach a proper calculation time. Present softwares 
and mesh generators allow to set up very complicated geometries 
(also polygonal), but it is a good point to start with simple mesh 
types and domain shapes (e.g. quadrangular, radial). A mesh 
refinement can be useful only for the area interested by mass/ 
energy transfers, like wells, natural recharge, atmospheric 
geothermal manifestation. Anyway the problem of calculation 
time saving must be considered, as it can occur for the operator 
to think about a massive refinement, on domains which are too 
big. Once the geometry and mesh have been set up, the equation 
parameters (and constitutive law) must then be considered. The 
constitutive law (e.g. Darcy Law) is typical for each phenomena, 


NUMERICAL 


SIMULATION 


RESOURCE 


= Type of reservoir 


= T, p, deliverability 


expected 


= Spatial extension 


ENERGY 
UTILIZATION 


Equations 
Model geometry 
Time/Space 
discretisation 
Mesh type 
Numerical methods 


NATURAL STATE 
(UNPERTURBED) 
SIMULATION 


UTILIZATION 
SCENARIOS 


I Q Power plant size 


QO) Mass flow rate 


SUSTAINABILITY 
ASSESSMENT 


Fig. 3. Information flow and prospect for the sustainable assessment through the use of the numerical simulation of the reservoir. 
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the conceptual model must arrange and properly describe the 
transport phenomena occurring in the geological system. Thermo- 
physical parameters should then be put into the simulation. It is a 
big problem, as the hypothesis of homogenous material is useful 
but it has some practical limitations. One single value of a 
significant parameter like porosity or permeability is assigned to 
a layer, or element (maybe with irregular shape). Often one has no 
possibility of assigning tensors or directional parameters (e.g. 
preferential directions for circulation), as they derive from very 
accurate measurements and interpretations. Calibration is a good 
tool, as it allows to address parameter values matching with real 
temperature and pressure data. Anyway it is possible only when 
reliable measured data are available. The main parameters to be 
assigned are essentially: 


© porosity (¢, defined as the ratio between the voids and the total 
volume considered) and effective porosity (taking into account all 
the “interconnected” voids, allowing then the fluid circulation); 

è permeability (k), being property of both the rock and the 
circulating fluid, giving an idea of the productivity of such a 
rock formation; 

© density (p, its distribution derives from accurate measurements 
and usually it is approximated with a single value for each rock 
formation); 

è thermal conductivity (of both rock and fluid), 

© specific thermal capacity (of both rock and fluid). 


The simulation is then run, and several well known algorithms 
can be used. The mathematical/engineering background is also a 
part of this very complex process. As the task is the sustainability 
of a project, a proper time scale should be identified, considering 
both economic lifetime (20 up to 50 years) and natural (or 
induced) renewability of the resource itself. Usually the steady 
(unperturbed) state is also simulated, to reach the conditions 
before the exploitation (in case of new fields). In this case the 
time scale is very large (10 up to 10” years). Large part of the 
models are used to simulate industrial scenarios of power/heat 
production and study the reservoir response during time. Different 
strategies can be adopted, in order to keep the utilization feasible 
and sustainable. It is very important to define the size of the 
installation (power plant, district heating) only after the simula- 
tion and complete resource assessment, this being a result of a 
sustainable exploitation strategy, not the starting point (Fig. 3). 
Oversizing or resource fast depletion (reservoir cooling, wrong 
reinjection) can be some consequences of incorrect sizing. 


3.2. Mathematical and physical background 


From a mathematical point of view, the studies about transport in 
reservoirs or aquifers were born in the field of mining engineering, oil 
and gas and applied thermodynamics. The task is to understand the 
response of a porous or fractured system when a hydraulic gradient is 
applied. Anyway the application of some experimental concepts to the 
geothermal fields is not a trivial task. The hypothesis and the real 
situation of the field should be compared very carefully. 

The basic calculations executed by the softwares for numerical 
simulation substantially deal with the resolution of the flow into 
porous and/or fractured media. Equations of conservation of the 
following properties are solved: mass, momentum, energy, con- 
centration of pollutant (or chemicals) dissolved. Some constitutive 
laws must also be implemented, dealing with the particular 
phenomenon involved. Geothermal reservoirs fluid flow can be 
generally studied according to the Darcy Law about porous media, 


the Darcy fluid velocity q is defined as (Saeid et al. [18]) 
k 
q= YP PENA) a) 


where k is the intrinsic permeability, being proportional to the 
hydraulic conductivity K according to the definition (in [18]) 


K = kg 2) 
H 


In which p is the density, g is the acceleration of gravity and p is the 
viscosity, p is the pressure and z is the vertical coordinate. One 
thing to be considered is that in case of anisotropy an hydraulic 
conductivity tensor K can be defined 


Kx Ky. Ke 
z 


R =| Kx Ky By (3) 
Kz Ky Ky 


This element is important in order to consider the complexity of 
the phenomena involved. Average scalar values of the hydraulic 
conductivity or permeability are often assigned, but the real behavior 
of rock formations can be very different from an averaged parameter. 
For example, porosity also changes with pressure and depth, it is 
important to understand which assumption must be done when 
assigning an average value to an extended rock body (having for 
example an unknown fractures field). The fractures are particularly 
difficult to reproduce in a numerical model, their accurate description 
is important if their spatial extension is to be compared to the rock 
formation size, otherwise average properties can be considered. In the 
most common reservoir simulation softwares, like for example 
TOUGHZ2 [19], it is also possible to define different models for the 
relative permeability (R;) as a function of the liquid saturation in the 
mixture (geothermal fluid) (Petrasim Manual [20]). 

Let us now consider a single phase (liquid) flow in a geothermal 
system. According to Rybach and Muffler [17] the following 
assumptions have to be considered: 


© rock matrix is homogeneous and isotropic, in particular refer- 
ring to porosity, permeability and thermal conductivity (sup- 
posed to be independent by temperature); 

è incompressible fluid; with density (p) and kinematic viscosity 
(v) dependent by temperature according to the laws: 


p =poll—a(T—To) — A(T — To)” (4) 
V=Voa(T) (5) 


with po, Vo, a, p and To being opportune constants, while o(T) is a 
function of T. Pressure work and dissipations due to viscosity can 
be neglected, internal energy of liquid (1) and solid matrix (r) being 


Ei = Civoi(T — To) (6) 


E; = Crvol(T — To) (7) 


with Cryo. and C;yo1 specific volumetric heat capacities of rock and 
liquid. Under these assumptions the balance equations of mass, 
momentum and energy can be written respectively as (according 
to [17]): 


vq =0 (8) 
(2) +ar—-t]1+40-T0] F +429 =0 9 
Po a K 
oT > 
[A — p)prcr + poc] zg trov d “VT =V-(kmVT) (10) 


km being the thermal conductivity of the mixture solid liquid. 


992 


A double phase system, in which the vapor phase is also 
considered, is described by similar equations [17]. The following 
assumptions about porosity have to be done in this case: 


© porosity ¢ only depends on local pressure of fluid; 
è liquid and vapor phases are in local thermal equilibrium. 


The equation of conservation (respectively) of mass, momen- 
tum (liquid and vapor) and energy can be then written as follows: 


o 
OOTY (ditrvdy)=0 (11) 
Rk 
d= ag (vp-m8 ) (12) 
Rk 
€ 
d=- 7 -(Vp—py £) (13) 
v 
0 
AlO- PPrEr + bE] + VET i+ PvE Fy +P di +P Gy) 
=V(kmVT)+(1 4 i+py dy) È (14) 


where E is the internal energy of the liquid-vapor mixture, K is 
the permeability tensor and R; is the relative permeability of the 
i-th phase (i=l, g liquid or vapor). Some fundamental differences 
between the flow in porous and coherent media than in fractured 
media must then be taken into account (|7,17]): 


(a) permeability induced by fracturing is generally higher than 
average formation permeability; 

(b) fracturing permeability is usually anisotropic (it depends on 
fracturing preferential direction); 

(c) the permeability due to fracturing is considerably more 
dependent on pressure and tension field in the rock with 
respect to rock matrix permeability. 


Fracture 


\ 
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In the recent paper by Zeng et al. [21] a discussion about the 
methods for fractured systems simulation is given. Geothermal 
fracture systems can be represented essentially in two ways: 


(1) discrete fracture network method (fracture orientation, spa- 
cing and mechanical properties); 

(2) equivalent continuous porous media method (the fracture system 
is seen as an equivalent continuous media, with similar meth- 
odologies as for double porosity and double permeability). 


In Fox et al. [12] an analytical and numerical model for 
fractured reservoirs (multi-fractured EGS) is presented and simu- 
lated. In the following Fig. 4 an example of simplified scheme of 
fluid circulation into the fractures is shown (b being the fracture 
size and x, the fractures mutual spacing). There is also the 
possibility of coupling thermal-flow problems with chemical or 
geomechanics problems into the same numerical simulation. It is 
the case of some models presented in literature, like for example 
the case of a fully-coupled flow and geomechanical model by Hu 
et al. [22]. 


3.3. Boundary and initial conditions 


The structural-physical model has to be set with the boundary 
and initial conditions, and also with constraints that can be useful 
for the results stability. The thermal boundary conditions (BC) 
usually represent the heat geothermal source entering the reser- 
voir: heat flow from the bottom, fixed temperature values at 
bottom/top or intermediate layers and adiabatic/impermeable 
conditions. BCs usually also represent the natural manifestations, 
natural recharge, lateral or regional flows, wells withdrawal/ 
reinjection in the aquifers and hydraulic head both for the 
hydrological problem and for the heat transfer. The BC kind are 
similar when considering flow or heat transport. Fig. 5 provides 


Fig. 4. Conceptual scheme for multi-fracture vertical geothermal system (x; is the spacing between fractures, b is the fracture aperture). 
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Natural outflow Natural recharge of fluid 


Inflow of hot fluid ==) = 
22 | 


Fixed value of temperature on the top 


\ 


Fixed value of temperature on the bottom 


Fig. 5. Example of boundary conditions (both thermal and mass conditions) in a 3D 
numerical model grid. 


example of boundary conditions in a 3D grid. Mainly four types of 
BC are used: 


(1) First kind (or Dirichlet) condition - along a border or a boundary 
the hydraulic head or temperature are assigned. 

(2) Second kind (or Neumann) condition - along a border or a 
boundary the fluid or heat flux is assigned. 

(3) Third kind (or Cauchy) condition — transfer coefficients are used 
particularly for the hydraulic head. 

(4) Fourth kind condition - single well or singular point source, 
typically used for wells conditions (extraction or reinjection, 
according to the sign convention), implementing Dirac 6 
function. 


In particularly dry reservoirs, conduction is the prevalent 
mechanism of heat transfer. In hydrothermal aquifers and tradi- 
tional geothermal systems convection flow also contributes to the 
mass/heat transport phenomena. Thermophysical parameters 
database are also available in the most used codes. Anyway, as it 
is stated in this work, a characterization of the parameters should 
by site-dependent in order to obtain reliable and physically 
consistent results. 

The initial conditions are usually the thermal gradient, pressure 
distribution in the domain and hydraulic head levels of rivers or 
reservoirs. To hold the results range near a specific value, con- 
straints about max/min fluid rate, pressure or heat flow can be set. 
Fig. 6 summarizes the various steps described here: first the 
definition of geological elements and dimension of the reservoir, 
then spatial discretization and materials calibration, finally the 
definition of boundary conditions, initial conditions and definition 
of temporal domain. 


3.4. Numerical integration of the equation systems 


The integration of all the interdisciplinary inputs and proce- 
dures is the most challenging and crucial part of a modeling 
process. An important issue of this process deals with the quality 
of information and data flowing through the starting phase and 
the simulation itself. Numerical simulation must be treated and 
used as an iterative process, continuously changing and improving, 


Geology — Spatial extension — Rock materials 


3D features - Spatial discretisation — 
Materials calibration 


Boundary conditions 
Initial conditions 
Time step (discretisation) 


Cold fluid 
(recharge) 


Post-processing, 
Results interpretation 


Fig. 6. Overall graphic view of the numerical simulation process (the second 
vertical geological section is from [23]). 


as the information flow goes both ways (Ungemach et al. [14]). 
A good and reliable model is often an object of discussion and 
reinterpretation. There are different techniques for space discreti- 
zation: finite difference, finite volume methods, as well as finite 
element methods. Different numerical integration techniques are 
implemented in codes and softwares. TOUGH2 [19] for example, 
uses a finite difference (time fully implicit) discretization, while in 
Feflow [24] finite element techniques are used. Other softwares 
are used in literature and in industry. Many softwares have been 
developed and used by Research Institutes or Universities. In those 
softwares a lot of the most known resolution algorithms (from 
numerical analysis and calculus) are implemented. The mesh 
refinement is a fundamental instrument that can be adopted to 
improve the analysis and optimize the computational tasks (con- 
centrating for example the mesh number in the wells area or along 
the faults). Different techniques can be adopted for modeling of 
the faults, but the accuracy about the data in input (upflow/ 
downflow, fluid rate, permeability, thermal anomaly) for these 
structures has to be very high to achieve reliable results. 

All the most used softwares are multipurpose, involving the 
possibility of simulating different types of diffusion phenomena. 
Pre and post-processors are typically embedded in commercial 
softwares, so that graphical interface and elaboration of the data 
can be easily carried out. The various design tools (also for shallow 
systems) are mainly based on the line-source model and Eskilson's 
approach, and the duct ground storage model DST (Haberl and Do 
[25]). When coupling of more phenomena is implemented it can 
give a very useful contribution for the scaling phenomena limita- 
tion (concentration, inflow-outflow, and variation during with- 
drawal-reinjection must be known). 
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Fig. 7. Typical dotted breakthrough temperature profile 


3.5. Limitations and criticalities of numerical simulation of 
geothermal reservoirs 


Notwithstanding its strategic importance it is clear that numer- 
ical simulation of geothermal reservoirs has some important 
limitations (which are treated and discussed in this work): 


© it strongly depends on the reliability and accuracy of the data; 
e numerically stable models can be physically inconsistent. 


The first limit can be also expressed by the principle “trash in- 
trash out”. It must be clear which is the physical-numerical 
problem to recourse numerical simulation for. One should evaluate 
if the numerical simulation is the more appropriate instrument to 
face a specific problem. Usually a problem can be simplified in a 
proper way to be solved according to calculus or numerical 
analysis without using dedicated softwares and complicated 


: (a) single extraction well; (b) and (c) double extraction well. 


geometric domains. Reservoir model simulation has to be pursued 
only if it is the most suitable and appropriate way to elaborate a 
design strategy. “Lumped parameters” models can be very useful 
for some medium-low temperature fields, considering plain litho- 
logical layers. Sometimes, particularly for linear and simple pro- 
blems they can be satisfactory, in spite of more sophisticated 
elaborations. 

One possible risk is to start “asking too much” or “asking too 
bad response” to the numerical models, making them “too much” 
or “too less powerful”. For example, starting from the same 
geological features of a field, a model can give different results 
depending on the resolution of the spatial distribution of the data. 
The numerical analysis could be important in order to define the 
exploitation scenario (Figs. 7 and 8) or the reinjection strategy. 

Some typical problems due to incorrect initial characterization 
of the resource can also be faced with an appropriate model 
simulation, they can deal with: oversizing of the plant (leading to 
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excessive extraction), scaling (causing corrosion, productivity 
drop, diameter reduction) and wrong reinjection strategy (losses 
of fluid or cooling of the reservoir). 


4. Numerical simulation of geothermal reservoirs: a review of 
case studies referred to experimentally investigated 
geothermal fields 


In this part of the paper various examples of numerical simula- 
tions of geothermal reservoirs are analyzed and critically reviewed. 
The point of connection among the various reviewed cases is the fact 
that they are referred to very well known geothermal fields for which 
a lot of experimental data are also available. Each simulation has its 
own peculiarity, concerning with dimensional scale, enthalpy level of 
the reservoir, fluid rate extracted/reinjected, and software used. A 
summary of the meaningful data is also given. All the cases described 
are different for various reasons. First of all for the typology of 
geothermal field: from medium enthalpy water-dominant field to 
dry-steam dominant field. The differences between the models deal 
with simulation domains (size ranges from some km up to about 
100 km), scenarios simulated (unperturbed or exploitation) and 
software used. 

One concept has to be emphasized in this review: the strong 
dependence of the results of the numerical analysis on the quality 
of the inputs and the difficulty that would be afforded in realizing 
the models. First of all, the data and the geo-thermo-physical 
parameters necessary are not always available or measurable. 

In the next two sections some cases will be analyzed: the review is 
organized as follows: in the first part are a group of cases related to 
geothermal field for which a lot of experimental data, boundary 


T (deg C) 


126 


42.8 


15.0 


Fig. 8. Graphical results about a simplified single well extraction system (also fluid 
velocity vectors and iso-temperature vertical section are shown). 


Table 1 


conditions are available and it is possible to reproduce the results of 
the simulation (Table 1). Each case is discussed in detail in a single 
subsection. A second group concerns cases available in the literature 
and refers to less investigated (and simulated) geothermal systems. 
The different techniques of numerical resolution of the equations in 
the model are also considered here. 


4.1. Momotombo reservoir (Nicaragua) 


The first case analyzed is the geothermal field of Momotombo 
in Nicaragua for which a meaningful set of data is available from 
literature. This case study has been chosen to remark how 
important the characterization and assessment of the geothermal 
potential is to undertake a correct industrial exploitation of a 
reservoir. This analysis is based on the literature data available, 
mainly from the work of Porras et al. [26-28]. 

Overestimation of the geothermal potential (and progressive 
plant oversizing) brought to a gradual and progressive impover- 
ishment of the resource in terms of temperatures, pressures and 
geothermal fluid rates. Since the initial years of industrial interest 
(1983-1989) the production decreased and different problems 
arised. A three-dimensional, porous-media, numerical model of 
the system was developed and calibrated and it was utilized to 
study the response of the geothermal reservoir under different 
fluid production and injection scenarios. The initial hypotheses 
are: total reinjection of the fluid extracted; reinjection at 100 °C 
and pressure 5 bars; time-constant fluid rates. 

The domain considered in this case is 3.1 x 2.4 km? with a total 
depth of —3000 m, divided into 9 layers. The whole domain is 
built with 972 blocks (minimum dimensions 200 x 200 m2, max- 
imum dimensions 600 x 700 m°). The total number of materials 
considered is 18. The most permeable layer is 5 x 107" m?, in a 
depth interval between —150m and —450m. Four different 
scenarios of exploitation, dealing with different extraction/reinjec- 
tion strategies are proposed. In [26-28] the authors state that 
satisfactory agreements were obtained between measured and 
computed discharge enthalpies and flow rates, for most of the 
shallow wells. The model also qualitatively reproduced the pres- 
sure drawdown measured in selected wells. 


4.2. Ngatamariki geothermal field (New Zealand) 


The Ngatamariki geothermal field is located 17 km from Taupo 
(North Island, New Zealand). It is one of the several high enthalpy 
geothermal systems within the Taupo Volcanic Zone (they are 
more than 20). Exploration wells were first drilled in 1985-1986, 
and Mighty River Power then drilled further wells in 2008-09. The 


Summary about the various numerical model analyzed: data referred to experimentally investigated reservoir. 


Reservoir Software Domain size and blocks no. 

Momotombo TOUGH2 3.1 x 2.4 km?; 3 km depth ; 972 
(Porras et al., [26-28]) blocks 

Ngatamariki TOUGH2 10.5 x 11 km?; ~5 km depth; 24128 


(Burnell and Kissling [29]) 
Larderello 1 (Romagnoli et al. [30], TOUGH2 

Barelli et al. [31], Arias et al. [32]) 
Larderello 2 (Della Vedova [33]) 
Wairakei (O'Sullivan et al. [34,35]; 

Mannington et al. [36,37]) 


blocks 


~ 10000 blocks 
SHEMAT 
Non-commercial 


codes; TOUGH2 8055 blocks (2006 model) 


Groß Schönebeck FEFLOW ~4.8 x 5.5 km?; ~0.6 km depth; 
(Blöcher et al. [39]) 489591 prism elements; 254744 
nodes 
Mt. Amiata (Barelli et al. [40]) TOUGH2 1100 km?; > 5.7 km (total thickness) 
Dubti (Tendaho Rift) TOUGH2/ ~2.5 x3 km?; >2 km depth 
(Battistelli et al. [41,42]) EWASG 


70 x 70 km?; ~7.5 km depth; 


42 x 26 km?; 10 km (total thickness) 
30 x 30 km?; 3.4 km (total thickness); 250 °C (Wairakei-Te Mihi) 


Mean res. T [°C] Fluid flow rate 


240-340 ~357 kg/s; 13 production wells; 
8 reinjection wells (power units: 32 MW) 


80 (scenario A); 120 ~695 kg/s (2 scenarios) 


(scenario B) 
200-300 °C 1300 kg/s (average historical data) 
200-300 °C (Only natural state simulation) 


~ 1460 kg/s (historical average); > 200 
wells (57 years) (power units: > 200 MW) 
~21 kg/s doublet of wells from the same 
drilling area 


150 °C (125 °C after 
30 years simulation) 


150-220 °C (first shallow (Natural state simulation; interaction 
res.); 300-350 °C(deep res.) between reservoirs) 
245-290 °C 140 kg/s 
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company has planned the installation of overall 130 MW electric 
power output, from a condensing unit and an ORC unit integrated 
in a combined cycle configuration (to be first run in 2013). The 
data about this geothermal field are available from a group of 
industrial reports by Burnell and Kissling (2009), in particular the 
one about the reservoir model [29], in which a very complete 
conceptual model can be found. The geothermal reservoir is 
composed of 3 parts: a shallow aquifer (50-100m deep); an 
aquitard in the rock formation of Huka Falls (50-400 m deep), 
between the surface aquifer and the intermediate one; and an 
intermediate aquifer, between Huka Falls rock formation and the 
clay cap. The grid has globally 26 horizontal layers (of different 
thickness), subdivided into 928 elements, reaching a total number 
of 24128 blocks. The surface extension of the model is 132 km?. An 
inlet of fluid flow (70 kg/s) is considered, entering from the bottom 
of the model, at a specific enthalpy of 1450 kJ/kg, with an heat 
flow of 101.5 MW. From the basement a constant heat flow of 
11 MW is considered. As constraints, areas with known flow 
direction and impermeable formations are simulated using fixed 
state (Teo, Pgeo) conditions. The validation of the model is based on 
temperatures and pressure data measured in different wells and 
on the calibration of other parameters like porosity, permeability, 
upflow mass and enthalpy, hydraulic connections between the 
deep reservoir and the groundwater system. 


4.3. Larderello geothermal field (Italy), 2010 model 


Larderello field (Italy) is one of the most anciently known and 
studied geothermal areas of the world. This field has been widely 
drilled and developed, with an almost 100-year-old history in 
geothermal energy utilization for power purposes. Average fluid 
production in the Larderello field, after the most recent explora- 
tions and improvements is now about 3700 t/h. The exploration 
extended in the early 1950s to the near Travale field (10-15 km SE 
of Larderello), which has now increased its fluid production to an 
average value of 1000 t/h. For this reason, when talking about 
large scale model of this geothermal system, usually one can talk 
about Larderello-Travale geothermal field. 

A numerical model about the field has been recently realized 
and improved, increasing the dimensions of the geological domain 
considered, by Romagnoli et al. [30], Barelli et al. [31], Arias et al. 
[32]. The model domain extent is 4900 km? (70 km each side), 
with a total thickness of 7.5 km. The grid is made of 10000 cells 
and 16 vertical layers. The geological scheme refers to five main 
rock formations: clayey-shaley caprock (0-500 m), fractured car- 
bonate reservoir (500-1000 m), metamorphic reservoir (1000- 
5000 m) and granitic intrusion as heat source of the system. 
Sixteen rock materials are considered and an impermeability 
condition along the boundaries is assigned. Fixed state (time 
invariant) conditions of temperature are considered at the top 
(15 °C, atmospheric pressure) and at the bottom of the producing 
layer (350-400 °C). Natural manifestations and cold inflow from 
the shallow aquifers are the only interactions with the external 
environment. The simulation of natural state has been carried out 
(millions of years as temporal scale) and a simulation of the 
exploitation history of the field has also been modeled. The 
historical data of 700 wells have been represented with 20 
“virtual” wells. The conclusion of Romagnoli et al. [30] are that 
only few changes in the conditions of the natural system have 
been caused by industrial development of the area. 


4.4. Larderello geothermal field (Italy), 2008 model 


A different model of this field has been proposed by Della 
Vedova et al. [33]. This model deals only with the natural state of 
the geothermal system, without considering the industrial 


exploitation. A very large temporal scale is considered (8-12 
millions of years). The extent of the considered domain are 
42 x 26 km’, with a total thickness of 10 km. The total depth of 
the model is very high, to include the K-horizons and the data 
from fluid inclusions. K-horizons are considered to be the main 
reservoir bottom, corresponding to the 400°C isotherm; collo- 
cated between 3000 m in the Larderello zone and 107m in the 
Travale zone. The numerical simulator SHEMAT has been tested 
and validated with the geophysical data available for the upper 
crust. The mesh cell size is 1x 1x0.3km?. The upper surface 
boundary conditions are 20°C and 0.1 MPa, impermeable and 
adiabatic conditions are assigned at the lateral boundaries. The 
bottom boundary is assumed to be impermeable and at a fixed 
temperature of 400-600 °C. A sensitivity analysis about thermal 
parameters is also considered in [33]. The authors remark that a 
lot of simulation have been run to match a composite target 
function, due to the uncertainty about several input data (geome- 
try, rock data). The work considered is an example of deep field 
simulation, oriented to the comprehension of the deep field 
phenomena more than to a sustainable exploitation approach. 


4.5. Wairakei geothermal field (New Zealand) 


Wairakei is another well-known geothermal field located in the 
Taupo Volcanic zone (New Zealand). The extent of the area is 
almost 25 km?. A complete report about the industrial history and 
the evolution of the numerical models can be found in O'Sullivan 
et al. [34] and Bixley et al. [35]. Fifty years of activity on this field 
have been reached in 2008. In the period 1958-1990 the power 
installed was increased up to 140 MW. After 1990 the power 
installed has approached 200 MW (with the plant of Poihipi, 
55 MW). The trend of the specific enthalpy of the geofluid in the 
field fits the evolution during the years [35]: rise of pressure due 
to reinjection activity; increase of temperature in the wells in the 
“Eastern Borefield” and seepage of cold water in the reservoir. 

Several models have been proposed and tested for Wairakei 
reservoir [34]. The models have been used to simulate different 
scenarios of field development. The common aspects and main 
characteristics of the conceptual schemes are: big fractures and faults 
(NE-SW), increasing permeability; two big upflow areas (260 °C at 
Wairakei, 300 °C at Tauhara); heat flux fixed values as BC, in the range 
300-600 MW; two regions at different pressures (high pressure at Te 
Mihi, low pressure at Southern Wairakei); natural recharge and 
rainwater infiltrations (1000 mm/year, 5% infiltration). A discussion 
about the calibration of the model is also reported in Mannington et al. 
[36,37]. The code iTOUGHZ2 [38], a TOUGH2 family software for inverse 
simulation, has been used to support the calibration process of the 
parameters and improve the match to the field data. Good matching 
results have been reached when the Tauhara field (strictly connected 
with Wairakei) has been introduced in the domain of the model. The 
perspective is to improve a 27886-blocks model (comprehending 
Wairakei-Tauhara—Rotokawa). 


4.6. Grofs Schönebeck reservoir (Germany) 


A case study about a simulation of an Enhanced Geothermal 
System (EGS) in a deep geothermal reservoir with hydraulic 
stimulation is discussed in Blécher et al. [39]. The geothermal 
research site of Groß Schönebeck is located 40 km north of Berlin 
(Germany). This case is presented as an example of well-based 
hydrothermal problem with a good match with the sustainable 
energy exploitation issue. A very good approach both to the 
numerical simulation and to the accuracy of the data is achieved. 

The software used for the simulations is FEFLOW [23]. The 
simulation deals with a doublet of wells. The reservoir is located 
between —3815 and —4247m below sea level, in the Lower 
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Permian of the Northeast German Basin. The hydraulic stimulation 
is a technology used to improve the productivity of a reservoir by 
inducing artificial fractures through high pressure fluid injection 
(water, gel-proppant). Six geological formations are considered, 
with different lithologies. Two sandstones formations (between 
—4000m and —4100m depth) are the most appropriate for 
geothermal power production. Natural fracture system is studied 
to analyze the conceptual scheme of circulation of the water, and 
to develop the induced fractures layout. The wells have a distance 
of 28 m at the surface, the injection well is vertical, while the 
production well is deviated to guarantee a sufficient distance of 
500 m between them within the reservoir. 

An interesting aspect remarked in [39] is the dependence of the 
hydraulic properties with temperature and pressure in the reser- 
voir (density, permeability, porosity, thermal conductivity) to be 
considered as part of the model. Also the total dissolved solids and 
chemical composition are considered in the whole model. The 
model has an extension of 4.809 x 5.448 km?, with a depth of 
almost 0.6 km. The grid is made by 489591 prismatic elements and 
254744 nodes, discretizing 27 spatial layers. The induced fractures 
are represented by 2D quadrilateral vertical elements. The produc- 
tion well has a fluid rate of 75 m?/h (~21 kg/s) at 150°C with a 
concentration of solids of 265 g/l. The injection well has the same 
fluid rate and concentration of solids but the temperature is 70 °C. 
The quasi-stationary conditions are reached after almost 4 x 10* 
years of simulation time (hydraulic head levels are matched). A 30 
years simulation of the exploitation with the above mentioned 
wells conditions gives a production temperature drop from 150 °C 
to 125.8 °C. 


4.7. Mt. Amiata geothermal system (Italy) 


In a recent paper of Barelli et al. (2010) the Tuscan geothermal 
system of Mount Amiata is considered [40]. The fields of Bagnore 
and Piancastagnaio have been explored and drilled (more than 100 
wells) for about 50 years. In this field two main reservoirs are 
present: the first one is in the carbonatic formations, between 500 
and 1000m deep, at average temperature of 150-220°C; the 
second reservoir is in the Paleozoic metamorphic basement at 
depths of 2500-4000 m, at temperatures of 300-350 °C. The Mt. 
Amiata Volcanic Complex has a total area of 80 km?. The peculiar- 
ity of the model is not only to analyze the exploitation of the 
reservoir but also to verify the possibilities of interaction between 
a phreatic aquifer (separated from the shallow aquifer by few 
hundred meters of impermeable formations) and the geothermal 
system. Moreover another particular aspect is that the two main 
reservoirs are characterized by gaseous caps (gas, vapor), in 
structures named “traps”. This is a peculiarity of Mt. Amiata field, 
occurring in the definition of the pressure distribution and fluid 
circulation model. The numerical model considered has been 
simulated with the software TOUGH2 [19]. The surface extension 
of the model is more than 1100 km?, with a total thickness 
between —4km and 1.738 km (a.s.l., Mt. Amiata top elevation). 
A time-constant heat flow (average 400 mW/m7?, with peaks of 
600-700 mW/m7) is the bottom boundary condition. The model is 
globally closed, referring to inflow-outflow of water. As remarked 
in [40], the outputs match with the field data. 


4.8. Dubti geothermal field (Tendaho Rift, Ethiopia) 


The model of the Dubti geothermal field review is based on a 
paper by Battistelli et al. [41]. The Tendaho Rift was identified as a 
promising geothermal area since an exploration project of the late 
1960s, and early 1970s. More recently (1990s) a drilling and wells- 
testing activity was carried out to explore and assess the geother- 
mal resource present in the area (central part of Northern Tendaho 


Rift). The drilling in the area of Dubti plantation confirmed the 
existence of a shallow geothermal reservoir. The temperature 
recorded in the drilled wells is in the interval 245-270 °C, while 
the temperature of the geofluid rising in a fault in the area is about 
290 °C (natural manifestations and fumaroles are present in the 
whole area). Lots of production/injection tests have been carried 
out in the area, and the paper cited is very detailed about these 
data and permeability distribution. The Dubti fault is considered to 
be ascending and upflow hypothesis is at the base of the 
conceptual model. Horizontal circulation (with eventually two- 
phase conditions), when crossing permeable layers, is also con- 
sidered. For the numerical model the software TOUGH2 [19] has 
been used, implementing the EWASG equation-of-state module 
developed for water, sodium chloride and CO, mixtures (Battistelli 
et al. [42]). The simple 3D model is extended about 7.5 km? with a 
total thickness of about 2 km. Natural meteoric recharges are 
considered (coming from the Ethiopian Plateau) together with 
horizontal and sub-vertical flows. The assumed initial temperature 
of the hot upflow is 290 °C. Model results and further drillings 
confirmed the presence of a hot fluid circulation zone at depths 
between 250 and 500 m in the Dubti area, also confirming the 
existence of more permeable zones along the Dubti fault. A 
possible development program is proposed in [41]: an extraction 
fluid rate of 140 kg/s from the shallow reservoir, for an expected 
power plant with maximum size 3-3.5 MW (back-pressure or 
ORC), serving the near region for about 50 years. 


5. Numerical simulation of geothermal reservoirs: an 
“extended” review of the other available cases 


In the present section a further review of several numerical 
simulations is reported. The cases discussed in the previous 
section and summarized in Table 1 deal with well known geother- 
mal fields and concern mainly with numerical modeling and data 
matching. The numerical simulations considered are referred to 
several geothermal fields object of analysis all around the world 
[21,23] and [43-70]. The additional cases listed in this section 
often deal with fields which are presented in the same paper 
together with the numerical simulation. Different softwares are 
used and various grid shapes and configurations are adopted. 

For some of the cases reviewed only a single paper is available, 
being the experimental data on the field not available, while in the 
most recent works a large amount of experimental and explora- 
tion data are available. Different geographic areas are involved, in 
order to perform a meaningful worldwide review. Also, different 
resource types are considered in the range of the medium-high 
enthalpy. The cases contained in this section have been reviewed 
in the Doctoral Thesis of one of the authors (Vaccaro [6]). Table 2 
illustrates the main characteristics of the geothermal field, namely 
the type of resource, fluid production, plant productivity (only 
estimated in some cases), reservoir extent and production depth. 
Table 3 contains the main indications about the numerical models, 
the objectives of the simulation, the geometry and extent of the 
model and grid and the boundary conditions. Some common 
aspects can be found in the boundary conditions and calibration 
process (due to the uncertainty about permeability). 


6. Discussion about the numerical models reviewed and 
perspective for utilization in the energy production perspective 


As it is evident from the previous review analysis, the numer- 
ical simulation of a geothermal reservoir is a well known field in 
the literature. In Sections 4 and 5 a summary of the cases is 
reported, they are all different for various reasons. First of all for 
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Table 2 
“Extended review”: characteristics of the geothermal fields. 


A. Franco, M. Vaccaro / Renewable and Sustainable Energy Reviews 30 (2014) 987-1002 


Field name and References Resource Production Plant: Size Extent Heat Max. 
(kg/s) and Type (MW) (km?) flux prod. 
(mW/ Depth 
m°) (m) 
Balcova-Narlidere, Turkey; (Gok et al. Water dominated, 140 °C 20.6° 0 z2 1100 
[43]) (annual 
avg.) 
Cerro Prieto, Mexico (Butler et al. [44]) Water dominated, double phase reservoir, 350 °C 620° 3000 
Den Haag, Netherlands (Mottaghy et al. Water dominated, single phase reservoir, 70 °C ~ 42° 5 (thermal, estinated) 547 60-70 2200 
[45]) (63) 
Desert Peak, USA (Zeng et al. [21]) EGS system, ~ 210 °C ~ 42 2-5 (EGS-derived 4.2 1200 
binary plant, under 
development) 
German Basin (Vogt et al. [23]) Water dominated reservoir, sandstone reservoir, 87 °C x42 District heating 25 50-78 2000 
Hammam Faraun, Egypt (Zaher et al. Water dominated reservoir, 70 °C Plant for desalinization 77 120 
[46] 
Hammam Musa, Gulf of Suez, Egypt Water dominated reservoir, 70 °C 13 (estimated potential, 45 80 
(Zaher et al. [47]) volumetric method) 
Hatchobaru, Japan (Yahara and Tokita Water dominated, double phase reservoir, 240-300°C  ~700 110 16.5 2300 
[48] (= 500 
reinjected) 
Heber, USA (Boardmann et al. [49]) Water dominated, double phase reservoir, 200 °C 120° (52 effective) 1830 
Hengill Area, Iceland (Gunnarsonn et al. Water dominated, double phase reservoir, 300 °C 330 210 2000 
[50] 
Kamojang, Indonesia (Suryadarma et al. Vapor dominated reservoir, 230-245 °C 300-424 200 14-21 ~ 1400 
[51]) 
Kizildere, Turkey (Yeltekin et al. [52] and Water dominated, double phase reservoir, 200 °C 250 (10 20.4 1240 
Ozkaya, [53]) effective) 
Las Tres Virgenes, Mexico (Guerrero- Water dominated, double phase reservoir, 250-275 °C 10-14 16.9- 800- 
Martinez and Verma [54]) 20.3 1000 
Los Azufres, Mexico (Maldonado et al. Water dominated, double phase reservoir, 240-280 °C, ~475 188 ~ 20 1500 
[55]) vapor dominated after exploitation 
Mt. Apo - Mindanao, Philippines Water dominated, double phase reservoir, 300 °C 104 2500 
(Esberto and Sarmiento [56], 
Emoricha et al. [57]) 
Pauzhetsky-Kamchatka, Russia Water dominated, double phase reservoir, 200 °C 250 (38 to 68 x5 63 1000 
(Kiryukhin et al. [58,59]) reinj.) 
Ogiri, Japan (Kunamoto et al. [60], 1st reservoir (400-200 m asl): water dominated, 50- x 320 (243 30 432 1800 
Itoi et al. [61]) 130 °C. 2nd reservoir (0 m asl): water dominated, double reinj.) 
phase, 230 °C 
Olkaria, Kenya (Ofwona, [62]) Water dominated, double phase reservoir, 200-300 °C 121 ~ 80 1600 
Onikobe, Japan (Nakanishi and Iwai 1st reservoir, water dominated, two-phase (—200 to 12.5 z1 175 1300 
[63]) 300 m a.s.l.), 200 °C. 2nd reservoir: water dominated 
(—1200 to —700 m asl), 250 °C 
Poihipi - Wairakei, New Zealand Vapor dominated 55.5 (19.5 55 (25 effective) 
(Zarrouk et al. [64]) reinj.) 
Sabalan, Iran (Bromley et al. [65], Water dominated, double phase reservoir, 140-250 °C 50+50 (planned) x100 ~200 zæ 3200 
Strelbitskaya and Radmehr [66], 
Noorollahi et al. [67,68]) 
Sumikawa, Japan (Pritchett et al. [69]) | Water dominated, double phase reservoir, 200 °C 0 400 2486 
Tanggu District (Guantao field), Tianjin, Water dominated, single phase reservoir, 60-75 °C 11-25 533 x 2000 


China (Lei and Zhu [70]) 


* The value is referred to a district heating system using the geothermal resource. 


P In 2000 an upgrade up to 720 MW was planned, a further upgrade to 820 MW was planned for 2012. 


€ Updated at 1996 [49]. 


the kind of geothermal field: from medium enthalpy water- 
dominant field to dry-steam dominant field. The differences 
between the models deal with simulation domains (size and 
features), scenarios simulated (unperturbed or exploitation) and 
software used. 

One concept has to be emphasized from this review: the strong 
dependence of the results of the numerical analysis on the quality 
of the inputs and the difficulty that would be afforded in realizing 
the models. First of all, the data and the geo-thermo-physical 
parameters necessary are not always available or measurable. 
Moreover defining the boundary conditions and initial conditions 
is not a trivial task. 

Initial conditions in the simulation are mostly: current thermal 
gradient measured; groundwater recharge flow in the reservoir 


and pressure distribution. The boundary conditions deal with both 
the hydro-geologic and the thermal problems, regarding essen- 
tially: hydraulic head; temperature distribution at the top/bottom 
of the geometrical domain; heat flow; natural recharge of steam or 
water; gases diluted in the geofluid and impermeable and adia- 
batic conditions at lateral faces. All of these conditions can vary 
with time during the simulated interval. 

In each case matching between the forecast and the measured 
data is then fundamental. The measured data must then be the 
basis for enhancing and improving the model during the following 
steps of the project. For most of the listed models and fields it is 
very difficult to report here a matching between modeled reser- 
voir and real measured data. In many cases the authors of the 
simulation are different with respect to the utilization plant owner 


Table 3 


“Extended review”: numerical model simulations. 


Field name and 
References 


Balcova-Narlidere, 
Turkey [43] 


Cerro Prieto, 
Mexico [44] 


Den Haag, 
Netherlands 
(Mottaghy 
et al. [45]) 

Desert Peak, USA 
(Zeng et al. 
[21]) 


German Basin 
(Vogt et al. 
[23]) 

Hammam Faraun, 
Egypt (Zaher 
et al. [46]) 


Hammam Musa, 
Gulf of Suez, 
Egypt (Zaher 
et al. [47]) 

Hatchobaru, Japan 
(Yahara and 
Tokita [48]) 

Heber, USA [49] 


Hengill Area, 
Iceland [50] 


Kamojang, 
Indonesia 
(Suryadarma 
et al. [51]) 

Kizildere, Turkey 
[52] 


Kizildere, Turkey 
[53] 


Las Tres Vírgenes, 
Mexico 
(Guerrero- 
Martinez and 
Verma [54]) 

Los Azufres, 
Mexico [55] 


Mt. Apo - 
Mindanao, 
Philippines 
[56] 


Simulation 


Natural state 


History matching (5 y) 
Evolution (20 y, 3 cases) 


Natural state 
History matching 


Evolution (30 y, 3 cases) 


Natural state 
Evolution (50 y) 


Natural state 
Evolution (20 y) 


Natural state 
Evolution (20 y) 


Natural state 

(10° y+3 x 10*y for 
fracture state) 
Evolution (20 y) 
Natural state 

(10° y+3 x 104 y) 


Natural state 


Evolution (9 y, 3 cases) 


Natural state 


Natural state (10 y) 


History matching (20 y) 


Evolution (30 y) 
Natural state 
History matching 


Evolution (30 y) 


Natural state 


History matching (17 y) 
Evolution (12 y, 9 cases) 


Natural state 
History matching 


Evolution (10 y, 4 cases) 
Natural state (7 x 10° y) 


Emplacement of 
magmatic chambers 
Evolution (30 y) 


Natural state 


History matching (25 y) 
Evolution (30 y, 2 cases) 


Natural state 


History matching (2 y) 
Evolution (5 y, 2 cases) 


Software 


TOUGH2 


TETRAD 


SHEMAT 
ArgusOne™ 
(grid) 


TOUGH2 


SHEMAT-Suite 


HYDROTHERM 


HYDROTHERM 


TOUGH2 


TOUGH2 
iTOUGH2 


STARS SAPHIR 
V.2.30 


SUTRA' v. 
1284-2D 


TCHEMSYS 


TETRAD 


TETRAD 


Geometry 


3D: 3x 4km?; 
thick.: 1.53 km (30 
to — 1500 m asl) 
3D: 11 x 10 km? 
thick.: 3.6 km 


3D: 22.5 x 24.3 km? 
thick.: 5 km 


vertical 2D-3D; 
height: 400 m; 
length: 400 m; 
thick.: 500 m 

3D: 5 x 5 km? thick.: 
6 km 


3D: 77 km? thick.: 
2.5 km 


3D: 45 km?; thick.: 
3.6 km 


3D: 16.5 km?; thick.: 
2.5 km (1100 to 

— 1400 m asl) 

3D: 14 x 13 km?; 
hick.: 3 km 

3D: 100 x 100 km; 
hick.: 2.9 km; 
(400 m to -2500 m 
asl) 

3D: 49.5 km?; 
hick.: 3.6 km 


3D: 840 x 600 km?; 
hick.: 1-1.2 km 
(estimated) 

3D: 870 x 720 km? 


3D: 20 x 30 km?; 
thick.: 22 km (1900 
to — 20000 m asl) 


3D: 12 x 12 km?; 
thick.: 3.2 km 
(3200-0 m asl) 

3D: 6 x 10 km?; 
thick.: 2.75 km 
(1250 m to -1500 m 
asl) 


Mesh 


Irregular grid 13 layers 5194 


blocks* 


7 layers x 1296 cells; 9072 
blocks‘; range size: 0.25- 


2km 


2.43 x 10° nodes 9 geol. units 


104 cells (x) 1 cell (y) 104 cells 
(z)(3x3,3x5,5x5) 


x 1.4 x 10° cells (18 
sedimentary layers) 


16 layers x (43 x 31) cells 


(length 100-1000 m) 


11 layers x (28 x 26) cells 


(length 100-1000 m) 


9 layers (thick.: 100-400 m) 


7425 blocks 


8 layers x 201 cells (1608 


blocks)" 


9 layers x 996 cells (8964 


blocks)* 


15 layers 12480 blocks 


5 layers (8 x 12 cells) 


29 x 24 cells 


I: 6307200; blocks II: 
4939200 blocks 


9 layers x (18 x 26) cells 


(4446 blocks)* © 8 


6 layers x (11 x 17) cells (1122 


blocks)* " 


Conditions and parameters 


Deep nat. recharge: 40 kg/s; shallow recharge: 11 kg/s; 
total heat rate input: 33 MW 


Top: T=40 °C; Lateral BC: T assigned to second and penultimate layers, 
lateral heat flux, fluid losses (347.2 kg/s); Bottom: T assigned, nat. recharge 
347.2 kg/s (350 °C); @ decrease with depth (range 0.176-0.01), fractures 
have constant porosity kx assumed as a logarithmic function of 

Top: constant T (11 °C, mean annual) Bottom: constant heat flux (63 mW/ 
m7); 4 depends on T 


Reference case: reinjection 1.5 kg/s (60 °C); ¢=0.2%; k=5 x 10° m?; Mass 
flow rate 50-75 kg/s 


Top: 11 °C surface temperature; Bottom: regional specific heat flow 
75 + 10 mW/m? (6000 m depth) Adiabatic lateral boundaries 
Top: T, p constant assigned; Bottom: heat flux fixed (120 mW/m?); Fracture 


zone conditions 


Top: T, p constant assigned (26.7 °C, 1.013 bars) 


Referred to a trial-error calibration process and 
further papers cited in [48] 


Top: T=25 °C, p=0,1 MPa; Lateral BC: impermeable, adiabatic; Bottom: T 
and p assigned, nat. recharge 

Top: T=15 °C, p=0.1 MPa‘; Penultimate layer: nat. recharge 1 kg/s 

(1500 kJ/kg) Bottom: T and p assigned ( ~ 265 °C), fluid flows 


T and p constant at boundaries (not specified) 


Boundary temperature at the value given by the gradient 
(0.03 °C/km), surface at constant 25 °C 


Top: natural emissions; Bottom: nat. recharge 8.3 kg/s (water at 345 °C) 


Natural emissions Bottom: nat. recharge 5 kg/s (320 °C), 10 blocks 


Calibration 


k” natural recharge 


k, 4? natural recharge 


Heat flux, 4 


Sensitivity analysis to ¢, k, 4, mass flow 
rate, Tyej, Po, productivity index 


heat flow, temperatures (3D inversion); 
Stochastic approach to hydraulic 
properties assignment 

Thermophysical properties determined 
by previous studies 


Trial-error calibration process 


k? 


k and natural recharge (iTOUGH2) 


¢, k (SAPHIR) 


T, p 


Rock parameters” natural recharge 


k, natural recharge 
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Table 3 (continued ) 


Field name and 
References 


Mt. Apo - 
Mindanao, 
Philippines 
[57] 

Pauzhetsky- 
Kamchatka, 
Russia [58] 

Pauzhetsky- 
Kamchatka, 
Russia [59] 


Ogiri, Japan [60,61] 


Olkaria, Kenya [62] 


Onikobe, Japan 
[63] 


Poihipi - Wairakei, 
New Zealand 
[64,65] 


Sabalan, Iran [66- 
68] 


Sumikawa, Japan 
[69] 


Tanggu District 
(Guantao 
field), Tianjin, 
China (Lei and 
Zhu [70]) 


Simulation 


Natural state 
(95 x 10° y) 


History matching (13 y) 
Evolution (5 y, 2 cases) 


Natural state 


History matching (36 y) 
Evolution (20 y, 3 cases) 


Natural state 


History matching (36 y) 


History matching 


Natural state (104 y) 


Natural state 


History matching (21 y) 
Evolution (10 y, 1 case) 


Prod. Wells 


Natural state 
Evolution (3 cases) 


Natural state (3 x 10* y) 
Evolution (50 y, 1 case) 


Natural state 
History matching 


Evolution (5 y, 2 cases) 


Software 


TOUGH2 


TOUGH2 
A-Mesh 


TOUGH2 
iTOUGH2 
A-Mesh 
TOUGH2 
iTOUGH2 


TOUGH2 


RANGER STAR 


TOUGH2 
iTOUGH2 
AWTAS" 
WELLSIM 
TOUGH2 


STAR 


AUTOUGH2 


* Mesh refinement in the production wells zone. 
> Manual iterative calibration. 
€ Matrix and fracture are both simulated in the cells. 


4 Both quadrangular and polygonal shaped blocks (the last ones are used to simulate the fracture). 


© Meteoric water recharge is also considered. 
f SUTRA - Saturated-Unsaturated TRAnsport. 
2 First layer reproduces orography. 


n Mesh refined at the sea level layers. 


Geometry 


3D: 22 x 26 km’; 
hick.: 3.25 km 
(1250 m to -2000 m 
asl) 

D: polygonal '; 
hick.: 0.7 km (avg.) 


w 


w 


D: polygonal’; 
hick.: 0.85 km (100 
o —750 m asl) 


3D: 5.5 x 3.9 km?; 
hick.: 2.85 km (250 
o —2600 m) 


3D: polygonal 
(120 km?)™; thick.: 
2.55 km (2000 to 
—550 m asl) 

3D: 6.5 x 8 km’; 
thick.: 2.4 km (400 
to —2000 m asl) 


3D: 12 x 8 km?; 
thick.: 4.6 km 


3D: 3 x 5 km’; 
thick.: 2.8 km (1200 
to — 1600 m asl) 
3D: 25 x 21.3 km’; 
thick.: 2.1 km 


i 7 sides polygonal area: 2 x 0.5 x 3 x 2 x 3x 1.75 x 1.5 km’. 


' MINC model applied to the main fault of the geothermal system (MINC - Multiple INteracting Continua). 


m 7 sides polygonal area: 11 x 11 x 5 x 9 x 10 x 2.5 x 6 km”. 
n AWTAS -Automated Well-Test Analysis System. 
° Mesh refinement in the atmosphere contact blocks. 


P First six layers reproduce orography. 


Mesh 


19 layers x (31 x 47) cells 
(27683 blocks, only 16411 


active)* 


131 cells x layer Double ¢ 


model 


3 layers (424 blocks, only 294 
active) Double ¢ model 


7 layers x (23 x 11) cells (1771 


blocks)* ! cell size: 


0.1-3 km (thick.: 0.1-1.6 km) 


5 layers x 158 cells (790 


blocks) 


14 layers x (10 x 11) cells 


(1540 blocks)* © ! 


16 layers x 192 cells; 2688 
blocks (2565 active)" °° ° 


16 layers x (9 x 10) cells 


(1440 blocks)°, P 


2 layers x (780 blocks) (thick.: 
150-1100 m); Block size: 
from 250 x 250 m? to 


2000 x 2000 m? 


Conditions and parameters 


Top: 65 °C; 0.1 MPa; Bottom: nat. recharge 145 kg/s (320 °C); Lateral BC: 
impermeable, adiabatic, outflow from constant pressure cells 


Top: nat. emiss. 100 °C; Lateral BC: constant T, p; Bottom: heat flux 63 mW/ 
m°, nat. recharge” 204 kg/s (830-875 kJ/kg), (SE) 120 kg/s (900 kJ/kg) 


Top: atm. T, p; nat. emissions (100 °C); Lateral BC: impermeable, constant T, 
p; Bottom: heat flux 63 mW/m?, nat. recharge 224 kg/s 


Top: T=75 °C, p=0.0981 MPa; Lateral BC: impermeable, adiabatic; Bottom: 
heat flux 43.2 mW/m? (432 mW/m?, S), total heat in 19.5 MW (260 mW/m? 
avg.); nat. recharge:° 30 kg/s (240 °C), total inflow (31.4 MW); 55 kg/s 
(1062.7 kJ/kg, inflow 58.4 MW), production area. 

Top: atm. cond.°, natural emissions (vapor) 366 kg/s; Lateral BC: E-W 
impermeable, N p=45 bars, S p=25 bars (1075 m asl); Bottom: nat. 
recharge 1253 kg/s (6 blocks), 1600 kJ/kg (avg.); fluid loss 958 kg/s. 


Top: const. patm, nat. emiss.°; Lateral BC: N-E impermeable and adiabatic; 
S-O constant p; Bottom: heat flux 175 mW/m?, 

nat. recharge 10 kg/s (330 °C) production area 

Porous/homogeneous media double ¢ (matrix/fracture) fractional 
dimension model 


Top: heat loss to atmosphere; nat. manifestations (modeled as 4 artificial 
wells, 40-50 kg/s). Second to last layer (near to the base): natural recharge, 
8 kg/s at 130 °C. Lateral: impermeable and adiabatic Bottom: nat. recharge, 
90 kg/s @ 265 °C (1159 kJ/kg, total 104.31 MW); uniform heat flux: 

200 mW/m? (Central-Southern area), null in the North area. 

Top: atm. cond. (T=10 °C), nat. Vapor emissions Lateral BC: E-W-S 
impermeable, adiabatic, nat. recharge; Bottom: impermeable, heat flux 
400 mW/m? (tot 6 MW) Initial: T vertical 10-250 °C 

Water table assigned. Top layer T and p constant (inactive elements). 
Natural recharge and inflows are assigned. 


Calibration 


k, @ natural recharge 


k, p natural recharge rock expansion 
coefficient 


iTOUGH2: T, p, natural emissions, 


natural recharge (m, k) History matchng: 


production wells", monitoring wells (T, 
P), $, k, fractures 
k (iTOUGH2) 


k? natural recharge 


k? natural recharge 


k, ¢ (AWTAS, iTOUGH2) 


K natural recharge 


k? recharge geology 


P thermophysical parameters” 


OOOL 
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while in other cases the history matching can be limited to a quite 
reduced time. 

Let us consider two different cases, the Larderello geothermal 
field and the Sabalan geothermal field. In the first case (see 
Table 1, [30-33]) a lot of historical data are available, also from 
detailed plant production, the simulation described in [30-32] has 
been realized to study the changes in the large field after years of 
exploitation and drilling, but it has been done only at the end and 
not in a forecast mode. For this very famous field it was impossible 
to think about a strategy as the one proposed here, the numerical 
simulation being a modern practical solution (the first plant in 
Larderello has been installed in 1913). In the case of Sabalan (see 
Tables 2 and 3, [65-68]) the numerical model is done during an 
exploration step, so the industrial production target proposed has 
not been verified yet. This is to discuss about the possibility, for 
some of the fields and models illustrated, of matching and 
comparing simulated and measured data. 

The different techniques of numerical resolution of the equa- 
tions in the model are not treated as a problem here. But it is clear 
that it is acommon problem with all the fields in which numerical 
simulation is involved. Although the codes used for this purpose 
are evolving very quickly and the results can be very detailed and 
widely complete, these simulations present a remarkable grade of 
uncertainty. In the perspective of a more diffused industrial 
development of medium to low temperature geothermal fields, 
the numerical simulation can be a very useful instrument that 
must be connected with the strategic elaboration and environ- 
mental and economic sustainability of the design of a geothermal 
plant. Overall a relatively good agreement was obtained in the 
various cases between the measured and computed temperature 
profiles, but the matching of pressure profiles appeared to be more 
difficult. Even if the definition of the model and of the boundary 
conditions requires particular attention and experience in order to 
avoid wrong results, numerical simulation could be a good 
strategy for an integrated design of geothermal plants and for 
the prediction of the geothermal reservoir response and environ- 
mental impact of geothermal plants [71]. 


7. Conclusions 


Numerical simulation is a fundamental and strongly interacting 
instrument for plant design. Different approaches to the numerical 
simulation of geothermal reservoirs operation are considered here, 
with reference to some case studies of geothermal fields and 
ground. The perspectives of numerical simulation of geothermal 
reservoirs as support to the design and sizing of geothermal plants 
are outlined. Models simulation is a powerful decision-making 
tool: it can provide useful indication about optimal sizing and 
sustainable management of geothermal utilization systems. 

The behavior of geothermal reservoir and time variations of 
temperature, hydraulic head and pressure have to be estimated 
before the design of the plant. This is a general statement for each 
utilization in geothermal energy (steam plants, flash plants, dis- 
trict heating systems) but it is particularly meaningful in case of 
binary plants based on Organic Rankine Cycles (ORC), whose 
efficiencies and operations strongly vary according to small varia- 
tions (concerning mass flow rate, temperature and pressure) of the 
characteristics of the resource. 

Moreover the issues of scaling and of correct definition of 
reinjection strategy must be considered. The importance of connect- 
ing geological-geophysical and energy engineering background 
appears fundamental for the success of the exploitation of a geother- 
mal field in order to sustain the geothermal fluid production rate over 
the whole lifecycle of the plant (in general higher than 15 years). 


A review of different cases of numerical simulation from 
literature is considered and discussed. The state of the art methods 
and commercial softwares available for the simulation of geother- 
mal fields are analyzed through the review of several geothermal 
fields and numerical models (ranging from medium temperature 
geothermal field to the well known cases of Larderello and Mt. 
Amiata, characterized by dry steam geothermal resource). Rela- 
tively good agreement was obtained in the various cases between 
the measured and computed temperature profiles; simulation of 
pressure profiles appears to be more difficult. 

Even if defining the model and the boundary conditions 
requires particular attention and experience in order to avoid 
wrong results, numerical simulation could be a good strategy for 
an integrated design of geothermal plants and for the prediction of 
the geothermal reservoir response after a long time exploitation. 
Nevertheless, the following issues have to be considered in order 
to obtain reliable results. 


e The geothermal phenomena simulation is complex, an inter- 
disciplinary approach is therefore necessary. 

© Software and codes often represent solver tools for the physical 
equations used. 

© Case study experience and history matching is a fundamental 
background, and its study should be enhanced. 

e Thermophysical parameters, boundary conditions, and mesh- 
ing method strongly affect numerical simulation. Only a high 
accuracy level of the input data provides reliable results. 
According to the quality of information available, simplified 
models can be adopted. 


An “integrated” approach to the complexity of the geothermal 
phenomena is still lacking. Geothermal energy is a particular 
renewable source: its use is sustainable only under particular 
conditions, which must be known particularly by investors and 
market players. 
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